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operator, the algorithm is free of systematic error and scales linearly with projection time and inter¬ 
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I. INTRODUCTION 

Quantum Monte Carlo (QMC) methods are power¬ 
ful and versatile tools for studying quantum phases 
and phase transitions. Algorithmic development in past 
two decades including the nonlocal updates^^® and the 
continuous-time formulations®’^ have greatly boosted the 
power of QMC methods, even surpassing the hardware 
improvements following Moore’s law. Using the mod¬ 
ern QMC methods, the simulation of bosons and unfrus¬ 
trated spin models is considered as a solved problem. 
QMC simulations therefore can be used to test novel the¬ 
oretical scenarios®^^^ and to verify experimental realiza¬ 
tions.^® 

While efficient algorithms exist for the simulation of 
bosons and unfrustrated spin models, simulations 

of fermions are more challenging than bosons and quan¬ 
tum spins because of the infamous fermionic sign prob¬ 
lem. It causes exponential growth of computational 
effort as system size or inverse temperature increases. 
Even for systems without a sign problem, the phase di¬ 
agram of correlated fermions can be nontrivial to estab¬ 
lish,®®’®® not to mention to accurately determine the uni¬ 
versality class and associated critical exponents.®®’®® The 
main reason for this difficulty is the unfavorable swperlin- 
ear scaling with system size and/or inverse temperature 
of determinantal quantum Monte Carlo methods, which 
are the workhorse of correlated lattice fermion simula¬ 
tions. 

Determinantal QMC method sums a factorially large 
number of fermion exchange processes into a matrix de¬ 
terminant, thereby avoiding the fermion sign problems in 
certain cases. The first algorithm based on this idea is the 
Blankenbecler-Scalapino-Sugar (BSS) method.®® It maps 
an interacting fermionic system to free fermions in a spa¬ 
tially and temporally fluctuating external field and then 
performs Monte Carlo sampling of this field. Numerical 


instabilities of the original approach have been remedied 
in Refs. 23 and 24. The BSS algorithm has become the 
method of choice of many lattice fermion simulations due 
to its linear scaling in the inverse temperature (5. We re¬ 
fer to Refs. 25 and 26 for pedagogical reviews. 

Closely related is the Hirsch-Fye algorithm,®®" which 
is numerically more stable and is more broadly applica¬ 
ble because it is formulated using a (potentially time- 
dependent) action rather than a Hamiltonian. How¬ 
ever, its computational effort scales cuhically with the 
inverse temperature and the interaction strength there¬ 
fore is much less efficient than the BSS method for the 
cases where both methods are applicable. The Hirsch-Fye 
method thus has typically been used in the study of quan¬ 
tum impurity problems and as impurity solvers in the 
framework of dynamical mean field theory (DMFT),®® 
where time-dependent actions need to be simulated. 

Both the BSS and the Hirsch-Fye algorithm are based 
on a discretization of imaginary time, thus introducing 
a systematic time step error, called the Trotter error. 
Nearly twenty years ago it was realized that the time- 
discretization is not necessary for the simulation of lat¬ 
tice models.®’®" Besides increased accuracy due to the 
absence of a Trotter error, continuous imaginary time 
formulations often results in a more efficient and flexi¬ 
ble algorithm.® In Ref. 29 a first continuous-time QMC 
method for lattice fermions has been proposed. How¬ 
ever the scaling of this algorithm and numerical stabi¬ 
lization have not been discussed in this paper and we 
are not aware of any application of the algorithm. Fur¬ 
ther development on fermionic continuous-time QMC 
algorithms®® have focused on quantum impurity prob¬ 
lems: the continuous-time interaction expansion (CT- 
INT) algorithm®®, continuous-time hybridization expan¬ 
sion (CT-HYB) algorithm®® and the continuous-time 
auxiliary field (CT-AUX)®® algorithm. CT-INT and CT- 
AUX are based on weak-coupling expansion of the action 
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TABLE I. Comparison between various determinantal QMC methods for fermions. The ground state methods are extensions 
of the corresponding finite temperature methods. They have similar scalings when replace the inverse temperature /3 by the 
projection time 0. N denotes the number of correlated sites and V denotes the interaction strength. 

Lattice Models Impurity Models 


Method name 

BSS 

— 

LCT-INT 

LCT-AUX Hirsch-Fye CT-INT CT-AUX CT-HYB 

Finite temperature 

Ground state 

Ref. 22 

Ref. 23, 34, and 35 

Ref. 29 

Ref. 30 

This paper — 

Ref. 27 

Ref. 36 

Ref. 31 

Ref. 37 

Ref. 32 

Ref. 33 

Trotter error 

Yes 

No 

No 

No 

Yes 

No 

No 

No 

Auxiliary field 

Yes 

Yes 

No 

Yes 

Yes 

No 

Yes 

No 

Scaling 

PVN^ ^ 

b 

l3VN^ 

PVN^ 

ipVNf 

{PVNf 

{PVNf 



^ Although the number of operations does not explicitly depend on the interaction strength V, one needs to increase the number of time 
slices proportional to V to keep a constant coupling strength with the auxiliary field, i.e. to retain the same level of fluctuations. 

The scaling of this code is unclear since it is not discussed in Ref. 29 and important implementation details are missing. 


and share the same scaling as the Hirsch-Fye method.^® 
These methods have revolutionized the simulation of 
quantum impurity problems and DMFT calculations.^® 
However, for lattice models they remained suboptimal 
compared to the BSS method due to their cubic scal¬ 
ing in the inverse temperature. Very recently an efficient 
continuous-time algorithm has been developed by two of 
the authors that scales identically to the time-honored 
BSS method®^ and can be used both with an auxiliary 
field (LCT-AUX) and without (LCT-INT). The prefix 
“L” indicating both their linear scaling and their appli¬ 
cability to lattice models. In Table I we summarize some 
properties of these determinantal QMC methods. 

Finite-temperature determinantal QMC methods can 
be extended to projector formulations,where the 
ground state is obtained from imaginary time projection 
of a trial wave function. In addition to being more direct 
to address quantum phases at zero temperature, the pro¬ 
jector formalism often allows for additional optimizations 
such as symmetry and quantum number projections^'^’^^ 
and combinations with fixed-node ideas in the presence 
of a sign problem."*^ In the case of the BSS method, nu¬ 
merical stabilization also becomes easier in the ground 
state formulation.On the other hand, for projection 
methods it is crucial to achieve a linear scaling in the pro¬ 
jection time since the results are exact only in the limit 
of infinite projection time.^® The ground state variants 
of the Hirsch-Fye and the CT-INT methods®®’®^ exhibit 
cubic scaling and thus are not ideal for lattice model sim¬ 
ulations. 

In this paper we present details of the projection 
version of the LCT-INT method whose feasibility has 
already been mentioned in Ref. 30. This algorithm 
provides an efficient continuous-time projection QMC 
approach for ground state simulations of correlated 
fermions. It retains the linear scaling with projection 
time and matches the one of the widely applied pro¬ 
jector BSS method®®^^®’^®’®^’®® while completely elim¬ 


inating the time discretization error. Moreover, the 
continuous-time formulation has greater flexibility for 
measuring observables and can easily be combined with 
histogram reweighting^^’"^® and extensive ensemble simu- 
lation®^’^®’"^^ techniques. 

The organization of this paper is as follows, in Sec¬ 
tion II we introduce a model system of spinless fermions 
that we will use to explain the algorithm in Section HI. 
Section IV contains comparisons of the method with 
other numerical approaches and results on the quantum 
critical point of spinless fermions on a honeycomb lattice. 
We end with discussions of future prospects in Sec. V. 

II. MODEL 

To make the presentation of our algorithm more con¬ 
crete, we will consider the following spinless fermion 
model at half-filling: 

H = Ho + Hi, (1) 

Ho = -tY^ (4ci + c]ci) = clK.jCj, (2) 
(ij) 

= (3) 

where Ci is the fermion annihilation operator, t denotes 
the nearest-neighbor tunneling, V > 0 denotes the ex¬ 
tended Hubbard repulsive interaction, and we have intro¬ 
duced the matrix K to denote the single particle matrix 
elements. 

Quantum Monte Carlo studies of this model on a 
square lattice date back to the early days of the BSS 
method.However, these simulations suffer from the 
fermion sign problem because the Monte Carlo weight is 
a single determinant which is not guaranteed to be posi¬ 
tive in general. Recently it was discovered that the model 
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(1) is naturally free from the sign problem on bipartite 
lattices at half-filling in the CT-INT formulation,®^’®^ be¬ 
cause the Monte Carlo weight can be expressed as the 
determinant of a real skew-symmetric matrix. This de¬ 
terminant equals the square of the matrix Pfafhan and is 
thus nonnegative. A conventional auxiliary field decom¬ 
position, on the other hand, breaks this symmetry. It was 
shown that this model also allows sign problem free sim¬ 
ulation in the BSS formalism if one works in a Majorana 
fermion representation,®^’®® i.e. performs the auxiliary 
field decomposition not in the density channel but in the 
pairing channel. The idea applies not only to the BSS 
algorithm but can be generalized to the continuous-time 
QMC algorithm. 


On the honeycomb lattice, this model exhibits a quan¬ 
tum phase transition from a Dirac semimetal to a charge- 
density-wave (CDW) phase. The quantum critical point 
is unconventional because of the coupling of CDW or¬ 
der parameter to the low-energy Dirac fermions.®^’®® 
Simulations using CT-INT found a critical point at 
Vc/t = 1.356(1) with critical exponents rj = 0.302(7) 
and V = 0.80(3).®® Although CT-INT is free from the 
time-discretization error, its cubic scaling with inverse 
temperature /3 limited these simulations to inverse tem¬ 
peratures pt < 20. To access the quantum critical point 
from a finite temperature simulation /3 was scaled linearly 
with the linear extent of the system, assuming a dynami¬ 
cal critical exponent z = 1. In Sec. IV B we will, as a first 
application of the projector LCT-INT algorithm, use it 
to directly address the quantum critical point of model 
(1) at zero-temperature and check our previous findings. 


III. ALGORITHM 
A. General Description 

In a projector QMC calculation one obtains the ground 
state wave function using imaginary time projection of 
a trial wave function [Tt) and calculate ground state 
observables as®®’®® 


^ (^T|e-e^|^r) ’ ^ 

For any I'I't’) with non-vanishing overlap with the true 
ground state, Eq. (4) approaches the ground state expec¬ 
tation in large 0 limit. In this paper, we choose |4'7’) as 
a single Slater determinant. 


Np / N \ 

= n |0), (5) 

j=l \i=l / 

where N is number of sites, Np is the number of particles, 
and P is a N X Np rectangular matrix.®® 

Instead of breaking the projection operator into 
small time steps as done in discrete time algo¬ 
rithms, ^^^^®’^®’®®’®® the continuous-time QMC formalism 
writes the projection operator in an interaction represen¬ 
tation 


e-e^ = Texp 



dr 


( 6 ) 


After a Taylor expansion of the exponential and time 
ordering the terms®®" the denominator of Eq. (4) reads 


(vI-T|e-®®^|vI/T)=5](-l)'= 

A;=0 






In the CT-INT and the CT-AUX methods,®® ’®2’®^ one 
applies Wick’s theorem to the integrand of Eq. (7) and 
expresses it as a determinant of a matrix whose size is 
proportional to the expansion order k. The subsequent 
simulation modifies the matrix with 0(kP) operations per 
Monte Carlo step. Since it takes k Monte Carlo steps to 
generate an uncorrelated sample, these methods®®’®^’®®" 
scale cubically with the average expansion order (fc). As 
{k) increases linearly with inverse temperature /3 (or 0 
in the ground state projection scheme) and interaction 
strength,®®’®^’®®" this unfavorable cubic scaling limits the 
applicability of these methods for lattice models at low 
temperature and strong interactions. Here we instead use 
the LCT-INT algorithm®® to achieve linear scaling with 


respect to the average expansion order. The algorithm 
scales as 0VA®, similar to the BSS algorithm. 

To proceed, we first express the interaction term 
Hi through exponentials of bilinear fermion operators. 
Traditionally this is done via Hubbard-Stratonovich 
transformation, at the cost of introduction of auxiliary 
fields. ^®’®^ Here we adopt a simpler approach based on 
the operator identity fii = |(1 — e™'^p. The interaction 
term Eq. (3) can then be expressed as 

(iJ) 

This reformulation works for any density-density 







4 


interaction®® and is crucial to respect the symmetry of lation. Substituting Eq. (8) into Eq. (7), the integrand 
the model (1), ensuring a sign problem free QMC simu- is recognized as a sum of determinants 


k=0 

det 


T^\ K /*0 

E E ••• E / dnj dT2...J dTfcX 


( 9 ) 


where K is defined in Eq. (2). The vertex matrix X{i,j) 
is an NxN diagonal matrix defined for nearest neighbors 
i and j whose non-zero elements are 




— 1, I = i or I = j, 
1, otherwise. 


( 10 ) 


Its form follows immediately from Eq. (8). In the follow¬ 
ing we denote the matrix product from imaginary time r 
to r' > r as a propagator 

= e-^^'-^-'>^X{irn,jm) ■ .. X (ti, , 

( 11 ) 

where ti and are the imaginary time locations of the 
first and the last vertices in the time interval. The cor¬ 
responding vertex matrices are X{ii,ji) and X{ijn-,jm) 
respectively. If there is no vertex in the time interval, 
the propagator simply reads = e~^'^ -r)K^ 

We then expand the nominator of Eq. (4) similarly and 
write the expectation value in a form suitable for Monte 
Carlo sampling, 


,A\ _ Sc ^(l^)(0)c,e/2 



( 12 ) 


where the configuration C denotes a point in the sum¬ 
mation and integration domain of Eq. (9) and (.. .)mc 
denotes Monte Carlo averaging according to the configu¬ 
ration weights w{C). 

For k vertices a configuration 


C = {{Ti;iiJi),{T2]i2,j2),---,{n-,ik,jk)}- (13) 


consists of ordered times 0 < ti < T2 < ... < < 0 

and corresponding pairs of nearest neighbor sites (*i, Ji), 
(* 2 i 72 ) ■ ■ ■ (ikijk)- An example of a configuration with 
three vertices is shown in Fig. 1. 

Using Eq. (11) the weight of a configuration is ex¬ 
pressed as 

w{C) = det[plp(0, 0)P] dri... dr^. (14) 

Since this weight is, up to a constant factor, identical 
to the weight of such a configuration in a CT-INT cal- 
culation®^’®^ at zero temperature,®^ these methods have 
identical sign problems. 

The quantum mechanical average {0)c,t of an operator 
O inserted into a configuration C at imaginary time r: 

{0)c,r = 

can be evaluated using Wick’s theorem since O is sand¬ 
wiched between two Slater determinants. 


B. Monte Carlo Sampling 

In this section we first explain how to sample using the 
weights Eq. (14) and then discuss efficient ways to per¬ 
form update and measurement using equal-time Green’s 
function. 



FIG. 1. A configuration with k = 3 vertices. We divide the 
imaginary time axis into intervals of size A and sweep through 
them sequentially. In each interval we propose updates that 
either insert or to remove a vertex. Measurement is performed 
at the center of the imaginary time r = 0/2. 


1. General Procedure 


To sample configurations C according to the weight 
w{C) we use the Metropolis-Hastings algorithm.®®’®® 
Starting from a configuration C we propose to move 
to a new configuration C with an a-priori probability 
A{C —?> C'). The new configuration is then accepted with 
probability p{C —i C') = min{l,r(C — C')} , where the 
acceptance ratio r is 


r(C ^ C') 


wiC')A{C' C) 

wiC)AiC-^C') ■ 


(15) 
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To facilitate fast computation of the acceptance rate, 
we divide the imaginary time axis into intervals of size 
A, as shown in Fig. 1. We focus our updates on one in¬ 
terval at a time, proposing several times to either insert 
or to remove an existing vertex at time t for site indices 
(i, j). Sweeping through the intervals we achieve ergod- 
icity. While such sequential updates violate the detailed 
balance condition, a global balance condition is still re¬ 
stored as long as the updates within each interval satisfy 
local detailed balance.®^ 

Using short hand notations, 

L{t) = and R{t) = B{t,0)P, (16) 

the insertion or removal changes R to R^ = X{i,j)^^R. 
Note that for the model Eq. (1) studied here i?+ = R~ 
because X{i,j)~^ = X{i,j), but the general Monte Carlo 
scheme does not rely on this property. The acceptance 
ratios are 


det(Li?+) VNbA 
det{LR) ^ 4(n -|- 1) ’ 
det(Li?“) An 
det{LR) ^ NbVA’ 


(17) 

(18) 


where Nb is the number of interacting bonds of the sys¬ 
tem, A is the length of the time interval on which we 
propose updates and n is the number of vertices in this 
interval. While insertion and removal updates are suffi¬ 
cient to ensure ergodicity of the sampling, one can nev¬ 
ertheless implement additional updates to improve the 
sampling efficiency, such as change the site index of a 
vertex (see Appendix A). 

After a full sweep through all the intervals, we measure 
the expection values of observables close to the center 
t = 0/2. 


2. Fast Update Using Equal-time Green’s Funetion 

Crucial for the performance of the algorithm is a fast 
calculation of the acceptance ratios Eqs. (17-18). They 
can be efficiently computed from the equal time Green’s 
functions Gim{T) = {cic\^)c,tj which in matrix form 
reads^^’^® 

G(T)=I-i7(r)[L(r)i?(r)]-iL(T). (19) 

The determinant ratio in Eqs. (17-18) can be expressed 
using the Green’s function as^^’^® 


det(Li?*) 

det(Li?) 


- det {I + [A(z, j-)^i - I](I - G)} = - det ^ ^ 

AGijGjt 


( 20 ) 


The second equality follows from [A(i,j)^^ — = 

—26iidim — ‘26ijdjm- With appropriately chosen trial wave 
function, the equal-time Green’s function of our model 
(1) has an important symmetry property which we prove 
in Appendix B: 


imaginary time t', which can be done by the following 
similarity transformation (assuming t' > r) 

G{P) = B{P,r)G{r)[B{r',r)]-\ (23) 


Gjiir) = Sij - ri^r]jGij{T), ( 21 ) 

where iji = ±1 for site i belongs to the A(B) sublattice. 
Similar to the case of CT-INT,®®’®^ Eq. (21) implies Ga = 
i and Gij = Gji for nearest neighbors. With this we 
can further simplify the determinant ratios to AGijGji = 
AG^j. Since the remaining factors in Eqs. (17-18) are all 
positive, there is no sign problem in the simulation of the 
model (1).®®’®! 

If a proposed Monte Garlo move is accepted, we update 
the Green’s function to G^ = I — R^{LR^)~^L using 


Q± — _ Glj{Gim 5 im) ^ Gli(Gjm ^jm) ^ 22 ^ 


G,;, 


G 




Propagating the Green’s function using the tricks dis¬ 
cussed in Sec. IIIG, Eq. (23) is more efficient than calcu¬ 
lating G(r') from scratch using Eq. (19). 


3. Observables 

All expectation values {0)c,t can be related to the 
Green’s function using Wick’s theorem. For example, 
the density-density correlation function can be expressed 

as {nifljri)c,T - (1 G/;)(l GyTjTn) -t- i^Slm Gjyil)GljYl- 

Using Eq. (21), the density-density correlation functions 
is measured as 


For a proof see Appendix G. 

In the Monte Garlo updates we always keep track of 
the Green’s function at the imaginary time r for which 
we propose an update. For the next Monte Garlo move 
we need to propagate the Green’s function to a different 












6 


the staggered CDW structure factor as 


M 2 


l,m ' ^ ^ 




(25) 


and the kinetic energy and interaction energy as 

(26, 

(i?i) = -^E ■ (27) 

{l,m) \ V ^ / / MC 

Another useful observable is the average expansion or¬ 
der 

. (28) 

/ MC 



Since there is no translational invariance along the imag¬ 
inary time axis, (k) is not directly related to the interac¬ 
tion energy (ffi) as it is in the finite temperature case.^^ 
Nevertheless, Eq. (28) still suggests (k) ^ 0VN, i.e. the 
average number of vertices scales linearly with the pro¬ 
jection time, the interaction strength and the system size. 
The fact that we are dealing with k of N x N matrices 
compared to the single 2k x 2k matrix of the CT-INT 
case^^ allows LCT-INT to achieve an 0{QVN^) scaling, 
as we will discuss in the next section. 


C. Algorithm Optimization and Complexity 

Achieving the same scaling of 0(QVN^) as in the BSS 
algorithm requires a careful implementation, for which 
an optimal choices of single particle basis and splitting 
imaginary time into intervals is crucial. 


1. Optimal Single-Particle Basis 

The main computational effort in performing the 
Monte Carlo updates is the propagation of the Green’s 
function to a new imaginary-time using Eq. (23) . Imple¬ 
mented naively this involves dense matrix-matrix multi¬ 
plication and requires 0{N^) operations, while the cost 
of the calculation of the determinant ratio in Eq. (20) 
is 0(1) and update of the Green’s function Eq. (22) is 
0(N'^). 

This unfavorable scaling can be circumvented by 
working in the eigenbasis of the noninteracting Hamil¬ 
tonian.^*’ In this way all the computations for MC 
steps Eqs. (23,20,22) can be performed with complexity 
0{N^). 

For this we have to use basis-transformed Green’s 
functions G = WGU, where U are the eigen¬ 
vectors of the single-particle Hamiltonian WKU = 


diag(Ei, £'2,..., i?Af). The basis change modifies the 
propagators to 

(29) 

{U^X(i,j)U)^^ = Sim - 2Ulu,m - 2Ulu,m- (30) 


In this basis the multiplication of G by either Eq. (29) 
or Eq. (30) requires only 0(N‘^) operations instead of 
0(N^).^° The disadvantage is that now the calcula¬ 
tion the determinant ratio Eq. (20) is slightly more ex¬ 
pensive. However since we only need one matrix ele¬ 
ment Gij = (UGW)ij which can be calculated using 
matrix-vector multiplication and vector inner-products, 
this 0(N‘^) overhead will not affect the overall scaling of 
the algorithm. Similarly, updating of the Green’s func¬ 
tion in the eigenbasis also keeps an 0(N^) scaling (see 
Appendix D). 

Working in the eigenbasis of the noninteracting Hamil¬ 
tonian does not increase the complexity of the measure¬ 
ments either. One can choose to measure single par¬ 
ticle observables in the eigenbasis and perform a ba¬ 
sis rotation afterwords. Alternatively one can rotate G 
back to G with 0{N^) operations for each measurement. 
Since measurements are performed only after a full sweep 
through all intervals, this does not affect the overall scal¬ 
ing of the algorithm. For many physical observables of 
interest we only need for neighboring sites (i,j) or 
for fixed site i because of translational invariance, which 
would reduce the required basis transformation to just 
0(N'^) operations. 


2. Optimal Interval Size 


Finally we show that by choosing the number of in¬ 
tervals M = 0/A proportional to the average number 
of vertices one can achieve an overall 0(QVN^) scal¬ 
ing in the algorithm. For each MG update, we need 
to propagate the Green’s function from some time t 
to another time r' in the same interval. This will on 
average pass through ^ (k) existing vertices, which 
is of order 0({k)/M). As we need 0{N‘^) operations 
to pass through each vertex and 0{N‘^) for calculat¬ 
ing the acceptance rate and for the actual update, one 
Monte Garlo step requires G(max{l, (k)/M}N‘^) opera¬ 
tions, where the max function accounts for the case of an 
empty interval. 

A sweep through all intervals and updating (k) vertices 
results in an overall number of G(max{l, {k)/M}N'^{k)) 
operations. By choosing the number of intervals M ^ 
(k) we can achieve an optimal scaling 0{{k)N'^) = 
0(QVN^). This should be compared to the scaling of 
other continuous-time methods which scale as (fc)^ ^ 

03y3jy3 31,32,37,38 
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D. Numerical Stabilization 

As in the BSS algorithm the multiplication of the 
Green’s function with the propagator for large 

imaginary-time suffers from numerical instabilities be¬ 
cause the matrix multiplication mixes large and small 
scales in the propagator. We stabilize the calculation fol¬ 
lowing a similar approach as used for the BSS method. 
The following discussion largely follows Refs. 24 and 26 
with the difference that our stabilization is done in con¬ 
tinuous time and in the eigenbasis of the single-particle 
Hamiltonian. 

To avoid accumulation of numerical errors we need to 
regularly recompute the Green’s function using G = I — 
WR{LR)~^LU, which requires us to have fast access to 
the matrices WR and LU in a numerical stable way. We 
thus divide the imaginary time axis into / intervals where 
within each interval the propagation is well-conditioned. 
The interval length is set by the inverse bandwidth and 
is independent of the total projection time 0. These 
intervals are different from the (shorter) intervals used 
for MG updates discussed above. 

At the interval boundaries (corresponding to imagi¬ 
nary times = 0,..., = 0) we store / -I- 1 ma¬ 

trices S^. Depending on the current imaginary time r 
of the Monte Garlo sweep, they hold the matrix product 
either to the right or to the left, 

(WRir*^), ifT>T^ 

(31) 

\L{t^)U, otherwise. 

On the rightmost and leftmost boundaries 5^^° = WP 
and = P^C/. The matrix 5'^ is updated whenever 

we cross the interval boundary in the sweep along the 
imaginary time axis. During a sweep from r = 0 to 0, we 
multiply the propagator = WB{t^ 

with to update 5^ = B{t\t^-'^)S^-\ In the 

backward sweep from r = 0 to 0, we update to 

We still need to stabilize the calculation of d 

themselves. Performing a singular-value-decomposition 
(SVD) 

U^R = UrDrVr, (32) 

LU = UlDlVl, (33) 

the different scales only appear in the diagonal matri¬ 
ces of singular values Dr and Dr. Since G = I — 
Ur{VlUr)~^Vl only depends on the well-conditioned 
matrices Ur and Vr, it is sufficient to keep track of them 
instead of the full matrix products. Therefore before up¬ 
dating we can perform an SVD on the matrix product 
P(r^ or B{t ^+^and only store Ur or 

Vl. 

Using these stored matrices - d ga,n easily re¬ 
compute the Green’s function at any imaginary time. To 
compute G(t) for > r > we can use the ma¬ 

trices 5'^+^ and and calculate Ur = B{t,t^)S^ and 


Vr = S'^+^P(r^+^,r). The Green’s function is then re¬ 
computed by G = I - Ur{VrUr)~^Vr. 

In the simulation we monitor the difference of the sta¬ 
bilized G and the old one to dynamically adjust the fre¬ 
quencies of the SVD and the recomputation of G. It 
turns out both frequencies are mainly set by the inverse 
bandwidth and are independent of the system size or the 
total projection time. Typically we need to perform one 
of such stabilizations for a propagation time r ~ 1/t. 
Since each of these stabilization steps costs 0{N^) due 
to the SVD or matrix inverse and we need to perform 
G(0) of them per sweep, it ends up with a scaling of 
0{QN^), conforming with the overall scaling of the al¬ 
gorithm. 

E. Calculation of the Renyi Entanglement Entropy 

Quantum information based measures play an increas¬ 
ing role in the identification of quantum phases and phase 
transitions.®^^®® In particular. Refs. 66-68 devised mea¬ 
surements of the Renyi entanglement entropy in deter- 
minantal QMG simulations. 

Since in the present algorithm the many-body ground 
state wave function is expressed as a sum of Slater deter¬ 
minants, the derivations of Ref. 66 concerning the re¬ 
duced density matrix hold. In particular, the rank-2 
Renyi entanglement entropy S 2 = — lnTr(/3^) of a re¬ 
gion A can be calculated using, 

-s, Ec,c' w{C)w{C') det [GaG (4 + (I - Ga)(I - G),)] 
^ Ec C' 

(34) 

where C and C are configurations of two independent 
replicas and Ga, G'^^ are the corresponding Green’s func¬ 
tion restricted to the region A. 

Although the estimator Eq. (34) is easy to implement, 
we observe it shows large fluctuations at strong interac¬ 
tion.®^’®®’^® We leave a discussion of extended ensemble 
simulations®^’®®’^^ in the LGT-INT formalism for a future 
study. 

F. Direct Sampling of Derivatives 

An advantage of continuous-time algorithms over 
discrete-time algorithms is that the Monte Garlo weight 
Eq. (14) are homogeneous functions of the interaction 
strength V. This allows to directly sample the derivatives 
of any observable with respect to V using its covariance 
with the expansion order k, 

= (§) + V ^ 

Higher order derivatives can be sampled in a similar way. 
Derivatives are useful for discovering quantum phase 
transitions and locating critical points. 






FIG. 2. Density-density correlation function of a 32-site chain 
with periodic boundary condition. Solid lines are DMRG re¬ 
sults. 


In discrete time approaches calculation of such ob¬ 
servables either relies on the Hellmann-Feynman theo- 
rem/^’^^ which is limited to the first order derivative of 
the total energy or requires noisy numerical dif¬ 

ferentiation of Monte Carlo data.^® 


IV. RESULTS 

We finally present results obtained with our algorithm, 
starting with benchmarks that demonstrate the correct¬ 
ness before presenting new results regarding the quantum 
critical point of the model Eq. (1). 


A. Benchmarks 

For all of our results we use a projection time 0t = 40 
and use ground state of the noninteracting Hamiltonian 
Hq as the trial wave function. In case of degenerate non¬ 
interacting ground states, we take as trial wave function 
the ground state of a system with anti-periodic boundary 
condition in a;-direction and periodic boundary condition 
along the y-direction. 

We start by showing results for a periodic chain. Fig¬ 
ure 2 shows the density-density correlation function of a 
periodic chain compared with results from density ma¬ 
trix renormalization group (DMRG) calculations, where 
C(r) is averaged over all pairs in Eq. (24) with the same 
distances r. At moderate computational cost we can per¬ 
fectly reproduce the exact ground state quantities using 
projection LCT-INT (filled dots). 

Figure 3 compares the results obtained with our algo¬ 
rithm to exact diagonalization (solid lines) for an V = 18 
site honeycomb lattice. Our method correctly produces 


FIG. 3. Interaction energy (red dots), ground state energy 
(yellow squares) and GDW structure factor M 2 (blue trian¬ 
gles) of an A = 18 site honeycomb lattice shown in the inset. 
Solid lines are exact diagonalization results. 




FIG. 4. Ground state energy per site of a honeycomb lattice 
versus inverse system length 1/L (QMC) or inverse bond di¬ 
mension 1/D (iPEPS). The QMG results of periodic boundary 
conditions and anti-periodic boundary conditions approach to 
the thermodynamic limit value from different sides. 

ground state results for the total energy, interaction en¬ 
ergy as well as the staggered density structure factor. 

Finally, we compare with infinite projected entangled- 
pair states (iPEPS) results obtained for the honeycomb 
lattice.iPEPS is a variational method which works 
in thermodynamic limit, whose accuracy can be system- 
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l/L or 1/D 

FIG. 5. The QMC results for the CDW structure factor com¬ 
pared with the square of CDW order parameter calculated 
using iPEPS. For V/t = 1 the CDW order parameter van¬ 
ishes for bond dimensions D > 8 in iPEPS. For V/t = 1.4 
we used a linear fit in 1/D of the CDW order parameter to 
obtain an estimate in the infinite D limit (see Ref. 51 for more 
details). 


atically improved by increasing the bond dimension D. 
Figure 4 shows the ground state energy per site versus 
l/L together with iPEPS results versus 1/D. QMC re¬ 
sults for systems with periodic boundary conditions and 
those anti-periodic boundary condition along a:-direction 
approach the L —>■ oo limit from different sides, thus 
bracketing the ground state energy in the thermodynamic 
limit. Extrapolation of all data yields consistent results. 
Figure 5 shows the CDW structure factor M 2 versus l/L, 
which extrapolates to the square of the CDW order pa¬ 
rameter. iPEPS on the other hand can directly mea¬ 
sure the order parameter since the symmetry is spon¬ 
taneously broken for an infinite system. Extrapolation 
again yields consistent results and shows the system or¬ 
ders for V/t = 1.4 but not at V/t = 1.0. 


B. Fermionic Quantum Critical Point 

We finally apply the projector LCT-INT to study the 
quantum critical point of spinless t-V model on a hon¬ 
eycomb lattice, which we previously studied by CT-INT 
simulations.®^ Our calculations go beyond the previous 
results in two aspects. We can directly address the T = 0 
quantum critical point using the projection version of 
LCT-INT and we are able to reach larger system sizes 



V 





FIG. 6. (a) Scaled CDW structure factor of different system 
sizes cross at the transition point (b) Scaled CDW structures 
factor collapse on to a single curve when plotted against scaled 
interaction strength. 

up to L = 18. Since a detailed finite size scaling study is 
beyond the scope of this paper, we use the critical values 
obtained in Ref. 51 and check for consistency. The CDW 
structure factor should follow the scaling ansatz 

M2L^+^ =f({V - Va)L^I'^^ , (36) 

where we previously found z + r] = 1.302, jz = 0.8 and 
Vc/t = 1.356.®^ Figure 6(a) shows the scaled CDW struc¬ 
ture factor M 2 L^~^^ where all curves cross around I 4 
when using these critical exponents. Scaling of the x-axis 
using (V — Vc)L^/‘^ yields good data collapse, shown in 
Figure 6(b). We conclude that the new zero temperature 
results on larger system size are consistent with previous 
findings concerning critical point and critical exponents 
in Ref. 51. 


V. DISCUSSION 

In this paper we presented details of the ground-state 
version of the LCT-INT algorithm of Ref. 30. As a 
continuous-time QMC algorithm it eliminates the Trot¬ 
ter error due to time discretization of the BSS algorithm 
while still keeping the favorable linear scaling with pro¬ 
jection time and interacting strength. It is therefore well 
suited for simulations of the ground state of strongly cor¬ 
related lattice fermions. 

Although the LCT-INT algorithms®® and the projec¬ 
tion version described here share operational similari- 
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ties with the BSS algorithm, 24 , 26,34 g^j-g ji^por- 

tant differences. In the BSS formalism, one breaks the 
projection operator into small discrete time steps 

and performs Trotter-Suzuki decomposition for each time 
step, which leads to a systematic time-discretization er¬ 
ror. The BSS algorithm then decouples the interaction 
terms using auxiliary fields. A typical update scheme is 
to sweep through these time slices^^’^® and flip the auxil¬ 
iary fields, similar to our scheme of sweeping through the 
intervals. However, the time slices of the BSS algorithm 
are fixed in time and their number is proportional to the 
projection time. Each time slice contains 0{N) auxil¬ 
iary fields, therefore even with a brute force propagation 
of the Green’s function on the site basis (Eq. (23)) one 
can achieve 0{N^) scaling. While in our case the number 
and positions of vertices are allowed to fluctuate so we 
need to propagate in the eigenbasis (Eqs. (29-30)) and 
use M ~ (fc) intervals such that on average each inter¬ 
val contains a single vertex to achieve a similar 0{N^) 
scaling. 


Formally the Monte Carlo weight in Eq. (7) is similar 
to the local weight of the CT-HYB method.In par¬ 
ticular the matrix version of CT-HYB^® also evaluates 
the Monte Carlo weight in the eigenbasis of a propaga¬ 
tor. However, our case is simpler because is a 

single particle propagator and the Monte Carlo weight 
simplifies to a determinant instead of than a trace in 
the many-body Hilbert space in CT-HYB. The present 
method can still benefit from algorithmic developments 
for the CT-HYB method. In particular a Krylov ap¬ 
proach for imaginary time propagation^® may bring the 
cost of propagation of the Green’s function to 0{N‘^). 
Our “sweep through intervals” scheme is also similar to 
the sliding window approach of the CT-HYB algorithm.®® 
One may alternatively consider using trees or skip list 
data structure to store partial matrix products.®®’®^ 
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Appendix A: Site-Shift Update 


An optional additional update shifts a vertex between 
X{i,j) by moving the site j to another neighbor of the 
site i denoted by j'. This will change the vertex matrix 
to X{i,j'). It amounts to insert a vertex matrix A(j, j') 
at the same imaginary time without changing the per¬ 
turbation order. The acceptance ratio is 


^shift 


det{LX{j,j')R) 

det{LR) 




(Al) 


Since the site j and f belongs to the same sublattice, 
Gjj! = —Gj'j (see Eq. (21)) ensures the acceptance ra¬ 
tio is positive. The formula for updating the Green’s 
function after a shift move is identical to Eq. (22), with 
indices i,j replaced by j, jb 


Besides being free of the discretization error, the 
continuous-time QMC approach provides a direct means 
to compute quantities such as observable derivatives 
Eq. (35), that are harder to obtain in discrete time 
simulations. These in turn may be used to locate in¬ 
teresting phase transitions with an accuracy that can¬ 
not easily be reached by standard discrete-time algo¬ 
rithms. Furthermore, the simple interaction dependency 
of the Monte Carlo weight w{C) ~ allows straight¬ 
forward combination with the histogram reweighting^^’^® 
or the Wang-Landau sampling^^’^® techniques. Both ap¬ 
proaches can produce results in a continuous range of in¬ 
teraction strength by recording histograms over the per¬ 
turbation order k. Combined with Eq. (35), the method 
offers new exciting opportunities to bring the study of 
quantum criticality of correlated fermions to a new level, 
approaching to what was achieved in the simulations of 
the classical®® and quantum spin systems.®® 


Appendix B: Proof of Eq. (21) 

Equation (21) is easiest to prove in the finite temper¬ 
ature formalism. Suppose the trial wave function I'I't) 
is the ground state of a noninteracting trial Hamiltonian 
Kt- The equal time Green’s function can be formally 
written as 

G(t)= lim ri-f H(r,O)e-^^"'H(0,T)l”\ (Bl) 

l3—>-oo 

We introduce a diagonal matrix = rjidij and the bi¬ 
partite conditions implies DKD = —K. Together with 
— ^ihj) this shows that (H(ri,T 2 )“®)^ = 
DB{ti,T 2 )D. Similarly, assuming the trial Hamiltonian 
also satisfies DKtD = —Kt, one has D. 

Combing these facts it is then straightforward to show 
that 



11 


Substituting Eqs. (C2-C4) into Eq. (Cl) and using that 
, for i ^ 7 , a~^ = 0 and ab = —l/(G,-,Gi,), one arrives at 

[G(r)]^-I= lim [l + S(0,T)^e-^^-S(r,O)^] -I Eq. (22). 

/3—)-oo i V / 


= - lim [I + (B(r, 0)-i)^e^'"- (i?(0, r)"!)^] 

/3->-oo 

= - lim \l + DB{T,Q)e-^^'^B{Q, t)D] 

^—>•00 '■ ■' 


-1 


-1 


= -DG{t)D. 


(B21 


Appendix D: Monte Carlo Updates in the Eigenbasis 


This proves Eq. (21). 


Appendix C: Proof of Eq. (22) 


Once a Monte Carlo update is accepted, we update 
G = WGU using 


Notice that X{i,j)im = (5;m(l - 2(5 h - 26mj), G^ can 
be obtained from G using the Sherman-Morrison formula 
twice:^'*’^® 


where 


GL = G'lm + bG[^{6,m 
G'lm — Glm + aGli{5im 



2 

" 1-2(1-G,,)’ 

2 

" 1-2(1-G',)' 


(Cl) 

(C2) 

(C3) 

(C4) 


Gfm = Glm - 


(guA (ug-u) 

iUGW)., ^ 

(gW) (ug-u) 

\ / li \ / jm 

iUGWh 


(Dl) 


This update only involves matrix-vector multiplication 
and outer-product of column and row vectors, both re¬ 
quiring 0{N‘^) operations. 
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